#--------------------------------------------------------------------------
# Tool Name:  Spatial Adjustment
# Source Name: SpatialAdjustment.py
# Version: 1.0
# Author: Hu Wang(wanghugigi@gmail.com)
#
#--------------------------------------------------------------------------
#ArcGIS 9.2 + PythonWin 2.6.4 + numpy + scipy

#differences between ArcGIS 9.2 and ArcGIS 9.3
#1. polygon.centroid as a point in ArcGIS 9.3 while a unicode string in ArcGIS 9.2

#--------------------------------------------------------------------------
#Imports
#--------------------------------------------------------------------------
# -*- coding: utf-8 -*-
import HelperFunctions as HF
import numpy as np
from numpy import array, random
from scipy.cluster.vq import vq, kmeans, kmeans2, whiten, py_vq2
from scipy import linalg
from datetime import datetime
from time import strftime


def getCurTime():
    """
    get current time
    Return value of the date string format(%Y-%m-%d %H:%M:%S)
    """
    format='%Y-%m-%d %H:%M:%S'
    sdate = None
    cdate = datetime.now()
    try:
        sdate = cdate.strftime(format)
    except:
        raise ValueError
    return sdate

def get_inputs():
    """Given the list of inputs for the script when run, return an object
    containing the variables for the script"""
    try:
        '''
        sInputFC = HF.pGP.GetParameterAsText(0)                     # Input Feature Class
        sZField = HF.pGP.GetParameterAsText(1).upper()              # Attribute Field Name
        sOutputFC = HF.pGP.GetParameterAsText(2)                    # Output Feature Class
        sXY = HF.pGP.GetParameterAsText(3).upper()                  # add the XY cordinate or not
        iNumCluster = HF.pGP.GetParameterAsText(4)                  # Number of Clusters
        lsField = sZField.split(";")
        iIter = HF.pGP.GetParameterAsText(6)                        # Number of iteration
        sMethod = HF.pGP.GetParameterAsText(5)                      # method of Kmeans
        '''
        '''
        #This block can be used to test the code using provided datasets
        sInputFC = "F:\Clustering by Python\countryClusters.shp"
        sZField = "CASES;Population"
        sOutputFC = "F:\Clustering by Python\cobc.shp"
        sXY = "No"
        iNumCluster = "10"
        lsField = sZField.split(";")
        iIter = 10
        sMethod = "Kmeans"
        '''

        #This block can be used to test the code using provided datasets
        sInputFC = "G:\Hu\Daofu\Daofu_County_Cluster03.shp"
        #sZField = "ElevMean;Temp;Rain;Yield"
        sZField = "ElevMean;Temp;Rain;Yield;Cluster_S"
        sCluster = "CLUSTER_S"
        sOutputFC = "G:\Hu\Daofu_SA.shp"
        sXY = "No"
        iNumCluster = "7"
        lsField = sZField.split(";")
        iIter = 2000
        sMethod = "Kmeans"
        
            
        
    except: raise HF.ReportError (HF.msgParseErr)
    
    #Make a generic object, populate it with our variables, then return it
    try:
        obj = HF.Gen()
        obj.add('sInputFC',sInputFC)
        obj.add('sZField',lsField)
        obj.add('sOutputFC', sOutputFC)
        obj.add('sXY', sXY)
        obj.add('iNumCluster', iNumCluster)
        obj.add('iIter', iIter)
        obj.add('sMethod', sMethod)
        obj.add('sCluster',sCluster)
        

    except: raise HF.ReportError (HF.msgInputsErr)

    return obj

def build_value_lists():
    """Reads through feature table extracting geometry and field info.
       Returns this information in dictionary structures."""

    pRows = HF.pGP.SearchCursor(inputs.sInputFC)
    # pRows = pGP.SearchCursor(inputs.sOutputFC)
    pRow = pRows.Next()
    sKey = np.array([])
    lsBadRecs = []
    iBadRecCnt = 0
    iRec = 0 

    HF.pGP.AddMessage ("--------" + getCurTime() + "--------")
    HF.pGP.AddMessage (HF.sBuildingLists)

    while pRow <> None:
        iRec = iRec + 1
        try:
            dTemp = []
            dTemp.append(pRow.GetValue("FID"))
            for Field in inputs.sZField:
                dTemp.append(pRow.GetValue(Field))
            dTemp.append(pRow.GetValue(inputs.sCluster))
            sKey=np.append(sKey,dTemp)
        except:
            iBadRecCnt += 1
        pRow = pRows.Next()

    #print iBadRecCnt     
    # Report any problems encountered reading feature input data.
    if iBadRecCnt:
        sMessage = HF.msgReadErr % (iBadRecCnt, iRec)
        HF.pGP.AddWarning (sMessage)
        HF.pGP.AddWarning(`lsBadRecs`)

    sKey.shape = (-1,len(inputs.sZField)+2)    
    #print sKey    
    return sKey

def calculate_certain_average(data, iClusterID):
    # this function is to find the average of the certain cluster which ID is ClusterID
    # data includes only attribute and cluster id
    # data does not include the index(FID)
    iLen = data.shape #iLen[0]:Line; iLen[1]:Row
    i = 0   #Line

    ndTempData = np.array([])

    while i < iLen[0]:
        if data[i,-1] == iClusterID:
            ndTempData = np.append(ndTempData, data[i,0:-1])
        i = i + 1
    ndTempData.shape = (-1,iLen[1]-1)
    average = np.average(ndTempData, axis = 0)
    average = np.append(average,iClusterID)
    #average = [average[0],average[1],...,iCluster,]

    return average
    
def calculate_all_average(data, iCluster):
    # this function is to find the average of each average in data
    # data includes only attribute and cluster id
    # data does not include the index(FID)
    # iCluster is a ndarray
    HF.pGP.AddMessage ("--------" + getCurTime() + "--------")
    HF.pGP.AddMessage ("Calculating average...")
    iLen = iCluster.shape
    i = 0
    ndAverage = np.array([])
    while i < iLen[0]:
        ndAverage = np.append(ndAverage, calculate_certain_average(data, iCluster[i]))
        i = i + 1
    ndAverage.shape = (iLen[0],-1)
    #print ndAverage
    return ndAverage


def calculate_deviation(data, average):
    iLenData = data.shape   #iLen[0]:Line; iLen[1]:Row
    i = 0   #Line
    
    iLenAve = average.shape
    
    dist = np.array([])
    
    while i < iLenData[0]:
        iClusterID = int(data[i,-1])
        j = 0
        while j < iLenAve[0]:
            if int(average[j,-1]) == iClusterID:
                temp = np.linalg.norm(data[i,0:-1] - average[j,0:-1])
                dist = np.append(dist,temp)
                break
            j = j + 1
        i = i + 1
    return dist
    

def calculate_single_deviation(data, average, iIndex):
    iClusterID = int(data[-1])
    i = 0
    iLen = average.shape
    #flag = 0
    #print iClusterID
    #print data

    while i < iLen[0]:
        if int(average[i,-1]) == iClusterID:
            #print "a"
            if abs(data[iIndex] - average[i, iIndex])< fRange:
                flag = 0
            else:
                flag = 1
            #print data[iIndex], average[i, iIndex], (abs(data[iIndex] - average[i, iIndex]) - fRange), flag
            break
        i = i + 1
    #print flag
    return flag

def output_results(Cluster):
    """ This function writes results to the output feature class. """
    
    # Add results field.
    HF.pGP.AddMessage ("--------" + getCurTime() + "--------")
    HF.pGP.AddMessage ("Adding results fields...")
    addField = []
    
    for Field in inputs.sZField:
        if len(Field) > 8:
            Field = Field[0:8]
        #Field = Field + "_N"
        addField.append(Field + "_N")

    if inputs.sXY == "Yes":
        addField.append("X_cent_N")
        addField.append("Y_cent_N")
    
    for Field in addField:
        if not properties.dcFields.has_key(Field.upper()):
            HF.pGP.AddField(inputs.sOutputFC, Field, "FLOAT")

    outField = "CLUSTER"    
    if not properties.dcFields.has_key(outField.upper()):
        HF.pGP.AddField(inputs.sOutputFC, outField, "LONG")
    else:
        outField = "CLUSTER_1"
        if not properties.dcFields.has_key(outField.upper()):
            HF.pGP.AddField(inputs.sOutputFC, outField, "LONG")
        
    #addField.append("CLUSTER") 

    # Add results to output FC
    HF.pGP.AddMessage ("--------" + getCurTime() + "--------")
    HF.pGP.AddMessage (HF.sWritingResults)
    pRows = HF.pGP.UpdateCursor(inputs.sOutputFC)
    pRow = pRows.Next()
    iError = 0
    iCnt = 0
    fInterval = len(attri) / 5.0
    fMore = fInterval
    iComplete = 20
    
    while pRow <> None:
        iKey = pRow.GetValue("FID")

        #print iKey
        try:
            if Cluster[iKey]:  # make sure we have a non-Null result.
                #since the type of Cluster[iKey] is numpy.int32, it should convert to [int]
                for Field in addField:
                    iIndex = addField.index(Field)
                    temp = float(clusterAttri_Norm[iKey, iIndex])
                    pRow.SetValue(Field, temp)

            temp = int(Cluster[iKey])
            pRow.SetValue(outField, temp)
                    
            pRows.UpdateRow(pRow)
            iCnt = iCnt + 1
            if iCnt > fInterval: fInterval, iComplete = HF.check_progress(fInterval, fMore, iComplete)
        #except: pass
        except:
            iError = iError + 1
        
        pRow = pRows.Next()
        

    #print "++++++++++++++++++++++"
    print iError
    
        
    HF.pGP.AddMessage (HF.s100Percent)
    HF.pGP.AddMessage(" ")
    pRows = None           
    #return sField  





def spatialAdjustment():
    import csv
    spContiguity = np.array([])
    #read contiguity file
    fn = "countyContiguity.csv"    
    for record in csv.DictReader(file(fn), dialect="excel"):
        spContiguity=np.append(spContiguity, [record['id1'], record['id2']])
    spContiguity.shape=(-1,2)
    
    def ContiguityFID(Fid):  
        temp = np.array([],dtype=np.int)
        iLen = spContiguity.shape
        i = 0

        while i<iLen[0]:
            if int(spContiguity[i,0]) - int(Fid) == 0:
                if abs(int(spContiguity[i,1]) - Fid) > 0:
                    temp=np.append(temp, spContiguity[i,1])
                #temp=np.append(temp, clusterAttri[int(spContiguity[i,1]),-1])
            elif int(spContiguity[i,0]) > Fid:
                break
            i = i + 1

        #temp = np.unique(temp.astype(int))
        return temp

    def ContiguityClusters(contiguity_fid):#the type of contiguity_fid is npaddary
        temp = np.array([],dtype=np.int)
        iLen = contiguity_fid.shape
        i = 0

        while i<iLen[0]:
            temp=np.append(temp, clusterAttri[int(contiguity_fid[i]),-1])
            i = i + 1

        temp = np.unique(temp.astype(int))
        return temp    
    
    def Contain(item, npList):
        iLen = npList.shape
        i = 0
        flag = 0
        
        while i < iLen[0]:
            if npList[i] == item:
                flag = 1
            i= i+1
        return flag

    def SingleConSameCluster(FID):
        FID = int(FID)
        tempCluster = ContiguityFID(FID)
        iCluster = clusterAttri[FID,-1]
        ndCSCluster = np.array([])
        
        for item in tempCluster:
            if clusterAttri[item,-1] == iCluster:
                ndCSCluster=np.append(ndCSCluster, item)
               
        return ndCSCluster

      

    def AllConSameCluster(FID):
        ndCSCluster = SingleConSameCluster(FID)
      
        for item in ndCSCluster:
            temp = SingleConSameCluster(item)
            for iitem in temp:
                if Contain(iitem, ndCSCluster) == 0:
                    ndCSCluster = np.append(ndCSCluster, iitem)
            if ndCSCluster.shape[0] >10:
                break
               
        return ndCSCluster
    
    #pRows = HF.pGP.SearchCursor(inputs.sInputFC)
    # pRows = pGP.SearchCursor(inputs.sOutputFC)
    #pRow = pRows.Next()
    #lsBadRecs = []
    #iBadRecCnt = 0
    #iRec = 0
    '''
    abc = 992
    print ContiguityFID(abc)
    print SingleConSameCluster(abc)
    print AllConSameCluster(abc)

    '''
    iLen = clusterAttri.shape
    i = 0
       
    while i<iLen[0]:
        #iRec = iRec + 1
        #try:
        #iFID = pRow.GetValue("FID")
        conCluster = ContiguityClusters(ContiguityFID(i))
        #print iFID
        if Contain(clusterAttri[iFID,-1],conCluster)==1:     
            if conCluster.shape[0]>1:
                ndSameClusterID = AllConSameCluster(i)
                if ndSameClusterID.shape[0] <10:     
                    if conCluster.shape[0]==2:
                        if conCluster[0] == clusterAttri[i, -1]:
                            iiClusterID = conCluster[1]
                            #print iiClusterID
                        else:
                            iiClusterID= conCluster[0]
                        #for item in ndSameClusterID:
                            #clusterAttri[item,-1] = iiClusterID
                        #elif conCluster.shape[0]>2:
                            

        else: 
            if conCluster.shape[0]==1:
                clusterAttri[i,-1] = conCluster[0]
            #elif conCluster.shape[0]>1:
                    
        #except:
            #iBadRecCnt += 1
        i=i+1

    #print spContiguity  
    
def output_results():
    """ This function writes results to the output feature class. """
    HF.pGP.AddMessage ("--------" + getCurTime() + "--------")
    HF.pGP.AddMessage ("Outputting results...")
    outField = "SpCluster"
    
    # Add results field.
    if not properties.dcFields.has_key(outField.upper()):
        try:
            HF.pGP.AddField(inputs.sOutputFC, outField, "text")
        except:
            HF.pGP.GetMessages(2)
    else:
        outField = "SpCluster_1"
        if not properties.dcFields.has_key(outField.upper()):
            try:
                HF.pGP.AddField(inputs.sOutputFC, outField, "text")
            except:
                HF.pGP.GetMessages(2)
      
                
    # Add results to output FC
    HF.pGP.AddMessage (HF.sWritingResults)
    #sFieldList = properties.sFID + ";" + sField + ";" + cField
    pRows = HF.pGP.UpdateCursor(inputs.sOutputFC)
    pRow = pRows.Next()

    # add message of processes 
    iError = 0
    iCnt = 0
    fInterval = len(attri) / 5.0
    fMore = fInterval
    iComplete = 20
    iKey = 0

    #test = []    
    
    while pRow <> None:
        iKey = pRow.GetValue("FID")

        #print iKey
        try:
            pRow.SetValue(outField, clusterAttri[iKey,-1])
                    
            pRows.UpdateRow(pRow)
            iCnt = iCnt + 1
            if iCnt > fInterval: fInterval, iComplete = HF.check_progress(fInterval, fMore, iComplete)
        #except: pass
        except:
            iError = iError + 1
        
        pRow = pRows.Next()
        

    #print "++++++++++++++++++++++"
    #print iError

    HF.pGP.AddMessage (HF.s100Percent)
    HF.pGP.AddMessage(" ")
    pRows = None

    #print test    
    #return sField  


#--------------------------------------------------------------------------
#MAIN

if __name__ == "__main__":
    """This function performs the 1995 Getis and Ord Gi* statistic."""
    
    inputs = get_inputs()
       
    if HF.pGP.exists(inputs.sOutputFC):
        HF.pGP.delete(inputs.sOutputFC)

    #Copy the input feature class to the output feature class.
    try:
        HF.pGP.QualifiedFieldNames = 0
        HF.pGP.Copyfeatures(inputs.sInputFC, inputs.sOutputFC)
    except:
        sMessage = HF.msgOverwriteErr % (inputs.sOutputFC)
        raise HF.ReportError (sMessage)

    properties = HF.get_featureclass_properties(inputs.sOutputFC)  
    attri = build_value_lists()
    #HF.pGP.AddMessage ("--------build_value_lists finished--------")
    
    clusterAttri = attri[:,1:-1]
    
    cluster = attri[:,-1]
    cluster.shape = (-1,1)
    clusterAttri = np.append(clusterAttri, cluster, axis = 1)
    iCluster = np.unique(cluster.astype(int))
    iCluster = np.sort(iCluster)   #iCluster is the unique cluster id in the cluster column, which is also sorted by ASC

    average = calculate_all_average(clusterAttri, iCluster)
    

    

    spatialAdjustment()
    #HF.pGP.AddMessage ("--------spatial adjustment finished--------")
    
    output_results()







    



#import csv

#reader = csv.reader(file('countyContiguity.csv', 'rb'))
#for line in reader:
    #print line
    #print line.count

'''
boards = []

for record in csv.DictReader(file('countyContiguity.csv'), dialect="excel"):  
    boards.append([record['id1'], record['id2']] )
print boards      
'''